Journal of Mathematical Biology

The global potential of a chemical reaction network has many applications and is closely related to the stochastic detailed balance. However, many fundamental questions concerning stochastic detailed balance remain unresolved, such as whether it depends on the system volume and how to construct new systems that satisfy it. In this paper, we show that stochastic detailed balance may depend on the system volume. We therefore introduce four types of stochastic detailed balance according to their dependence on volume and rate constants, and systematically investigate the relationships among them. Our results distinguish detailed balance arising from particular choices of volume and parameters from that enforced by network structure, and identify conditions under which detailed balance at one volume extends to all volumes. We further obtain a class of networks satisfying stochastic detailed balance for every volume and every positive choice of rate constants, and construct new systems whose global potentials exhibit double-well structures.

AO_SCPLOWBSTRACTC_SCPLOWCellular automata and graph reaction-diffusion systems encode local spatial interactions in different mathematical forms. We develop a cochain-operator calculus for these two settings. Over a finite field Fq, every local rule on a finite neighborhood has a unique reduced polynomial representative. On an oriented line, the coboundary and endpoint maps recover the left and right shifts. Our main theorem shows that these operators, together with linear operations, constant cochains, and the degree-zero cup product, generate every finite-radius polynomial cellular automaton. Explicit formulas for Rules 30, 110, and 22 show how reflection-invariant linear coupling, directed transport, and nonlinear neighbor interactions enter the calculus. On a general graph, d*d is the unweighted combinatorial Laplacian and enters a graph reaction- diffusion recurrence. Over [R], the term - Dd*d with D [≥] 0 admits the usual diffusion interpretation; over Fq, the corresponding expression defines modular coupling without an intrinsic order. In the morphogenetic examples, we therefore distinguish pattern-generating dynamics from finite-state observation and use the Betti numbers of active induced subcomplexes to summarize observed patterns. This yields a common algebraic representation without identifying real-valued diffusion with finite-field dynamics.

Fractional differential equations are increasingly used to represent temporal memory in tumour immune and pharmacokinetic models, but changing the derivative order without a consistent time scale can make comparisons difficult to interpret. We study a dimensionless, commensurate Caputo system for tumour, immune, normal cell, and drug states. We establish positivity, boundedness, the tumour invasion threshold, and local stability conditions. We then use a fractional Adams Bashforth Moulton predictor corrector method with an explicit mesh convergence study. Two conclusions correct potentially misleading interpretations of the model. First, the equilibrium locations are independent of the fractional order: decreasing the order changes transient time scales but does not by itself alter the model's asymptotic states. Second, after a unit bolus the drug equation has the exact solution C({tau})=E(-K{tau}). Its maximum is therefore the initial concentration, not a delayed peak. For 0<<1, the algebraic tail makes the infinite-horizon area under the concentration curve diverge. In the illustrative coupled simulations, the apparent late tumour reduction at small is consequently driven by persistent modelled drug exposure and must not be read as evidence of superior therapy. All parameters are dimensionless and illustrative; no patient data are fitted. Source code and generated tables are included to make every numerical statement reproducible.

Non-exponential residence times give rise to history-dependent or memory effects in epidemic dynamics, yet their mathematical formulations remain fragmented. Here we develop a unified framework for a history-dependent SEIR model with general latent- and infectious-period distributions. We first establish the equivalence between the stage-age and renewal formulations and show that both admit a unified transition-flow representation. We then introduce the concept of a transition-flow multiplier, which provides a direct link between residence-time distributions and transition flows. When the multiplier is a weighted sum of shifted power functions, the model admits an exact representation in terms of Riemann-Liouville or tempered Riemann-Liouville fractional operators, providing a first-principles foundation for fractional epidemic models. Finally, using Euler-Lotka equations and stochastic-order arguments, we characterize how the latent- and infectious-period distributions shape epidemic growth and decay. A longer latent period slows both early growth and late decay, while, at fixed means, greater variability in the latent period accelerates early growth but slows late decay, whereas greater variability in the infectious period slows both.

Moment-closure methods provide a widely used approach for approximating low-order statistics of stochastic biochemical reaction networks, but the resulting closed equations may exhibit physically meaningless dynamics. Here, we systematically investigate the stability of second-order moment-closure approximations for a bursty gene expression model with nonlinear protein degradation. We consider the Gaussian, log-normal, and negative-binomial closures. By analyzing the resulting dynamical systems, we show that the Gaussian closure fails the proposed stability criterion for all parameter values, whereas both the log-normal and negative-binomial closures possess a unique physically accessible fixed point and are stable if and only if kB(1 + 2B) [&ge;] 1, where B is the mean burst size and k is the dimensionless protein synthesis rate. Moreover, this condition is automatically satisfied when the steady-state mean protein copy number is at least one, suggesting that it is mild in biologically relevant regimes. Numerical comparisons over a broad parameter range further show that the negative-binomial closure is the most accurate, while the Gaussian closure is the least accurate. We further develop higher-order negative-binomial moment-closure schemes and demonstrate that increasing the closure order substantially improves predictive accuracy. These results highlight the importance of considering dynamical stability, in addition to numerical accuracy, when assessing moment-closure approximations.

Much of the literature and intuition associated with mathematical epidemiology is driven by deterministic models, which are a reasonable assumption when the population size is large. Stochastic models, especially individual based models, are however considered vital when dealing with small population sizes, especially at times of invasion or extinction. The overwhelming majority of these models (both deterministic and stochastic) assume that the underlying parameters are fixed (or follow a regular seasonal pattern). Here, we consider an analytic framework for dealing with randomly varying parameters through the use of stochastic differential equations - thereby capturing the action of external noisy processes such as weather. In particular, we focus on when the transmission rate, {beta}, varies as the solution to a Cox-Ingersoll-Ross Model, such that {beta} is gamma distributed with autocorrelation. We consider the impact of this parameter variation on a stochastic version of the Susceptible-Infected-Recovered model, and for this 'double-stochastic' model show through simulation and analytical results that exogenous noise increases the impact of stochasticity, potentially leading to more early extinctions, wider variations in the number of cases at equilibrium, but that early growth rate can be faster or slower depending on the precise parameters.

The dynamics of gene regulatory networks are governed by intrinsic noise, stemming from the random nature of biochemical reactions, and by extrinsic noise, arising from fluctuations in cellular components and environmental conditions. Together, these sources can compromise the reliability of predictive computational models if not properly accounted for, and capturing both effects within a single framework remains a non-trivial task in computational biology. In this work, we propose an uncertainty quantification framework that addresses these two contributions jointly: intrinsic stochasticity is described through a partial integro-differential equation (PIDE) for the protein probability density function, whereas extrinsic noise is represented as parametric uncertainty in the kinetic parameters. The propagation of the uncertainty is carried out via an intrusive polynomial chaos expansion (PCE), in which the PCE coefficients are obtained from a stochastic Galerkin projection of the PIDE, yielding a coupled deterministic system that is solved with standard numerical methods. We illustrate the approach on a positive autoregulatory gene network with one and two uncertain kinetic parameters. The proposed approach accurately reproduces the mean, variance, and full protein probability density function, including the bimodal distributions, at a substantially lower computational cost.

We show that the p53-Mdm2 gene interaction can be understood as a form of antithetic feedback, but with mixed feedback loops. In particular, we show that this viewpoint sheds light on the emergence and robustness properties of p53's oscillatory behaviour, as observed experimentally.

The super-organismal view of ecosystems has largely been superseded by the individualistic view. One consequence of the dominance of the individualistic view is that many modern ecologists treat ecosystems as nothing more than vectors of relative abundances, ignoring the potentially important idea that an understanding of ecosystems should be based on dynamics rather than abundances or species identities. We develop a mathematical framework in which we compare dynamical properties of ecosystems with different sets of species, using ideas from functional analysis, metric spaces, and topology. We give two proof-of-principle applications of our framework to marine sessile communities and to a large database of ecosystem models. We show that under a set of biologically-motivated axioms designed to capture the properties of predator-prey systems, there is only one natural kind of ecosystem.

We propose a minimal mathematical framework to describe the go-or-grow dynamics of tumor cells comprising two phenotypically distinct populations. One population is migratory and undergoes constant diffusion, while the other proliferates in an oxygen-dependent manner. The local oxygen concentration governs transitions between these phenotypes. We then ask whether the time scale of cyclic hypoxia, where oxygen is spatially homogeneous but fluctuates between normoxic and hypoxic levels over time, alters the dynamics of the reduced generalist model in the fast phenotype-switching regime compared with the full go-or-grow model. We observe that the reduced generalist phenotype model fails to capture the go-or-grow dynamics under rapid oxygen fluctuations. Our results predict that the tumor growth rate increases in the generalist-phenotype model under rapid oxygen fluctuations compared with tumors of the proposed go-or-grow model. We then explore whether oxygen consumption, which leads to spatial heterogeneity in oxygen concentration, plays a role in distinguishing the dynamics of the reduced model from the proposed go-or-grow model. We demonstrate that the tumor in the proposed go-or-grow model exhibits a higher growth rate than the tumor in the reduced generalist-phenotype model under a slow oxygen consumption rate.

Microbial growth has often been analyzed under the assumption that microorganisms have evolved to optimize phenotypic characteristics of interest, such as growth rate and growth yield. This assumption has been useful, for example, for the genome-scale modeling of metabolism and the study of the allocation of cellular resources to physiological processes. In many experimental situations of interest, however, microorganisms are found to be suboptimal with respect to phenotypic characteristics that are thought to be favorable in that situation. Whereas a strong theoretical framework exists to mathematically relate cellular resource allocation strategies to Pareto optimality of microbial growth and other biological processes, much less is known about the consequences of Pareto suboptimality. We extend the framework to the latter case and show that a given Pareto suboptimal phenotype can be explained by a range of underlying resource allocation strategies, each corresponding to a different growth physiology and biomass composition. We test the predictions with the help of a coarse-grained model of microbial growth and published experimental data, which relate Pareto suboptimal rate-yield phenotypes of Escherichia coli and the microalga Tisochrysis lutea to the macromolecular composition of the cells (storage metabolite and total protein contents). Changing the focus from Pareto optimality to suboptimality provides interesting leads to exploring the diversity of growth strategies that can support a given phenotype. This change of perspective is of practical interest, because some of the growth physiologies within this range may be important for biotechnological applications

The study of tissue dynamics has been stimulated during the last decades thanks to the use of quantitative descriptions, with the development of several theoretical and computational frameworks, many of them revolving around the notion of reaction-diffusion systems, eventually with additional structure variables beyond time and space. The use of structure variables can accommodate phenotypic traits. In this work, we study a family of competition models, where a given population depends on a resource (e.g. oxygen) and several populations are competing for it. Our quantitative description incorporates phenotypic traits and heterogeneity at the level of cell cycle variations, which influence replication rates via oxygen consumption. This enables us to replicate the fitness of specific subpopulations to environmental conditions (e.g. oxygen shortage or external influences). Using numerical simulations, we show that such models display dynamical pattern formation in the form of coupled travelling wave profiles that expand or retreat at the same wave speed. The full theoretical analysis of such dynamics is quite involved; to circumvent this difficulty, we introduce a quasi-stationary approximation for the resource dynamics. We find that this approximation can reproduce the overall behaviour very accurately, with the additional benefit of allowing theoretical treatment of the reduced model. In this way, we provide estimates on the wave speed which are numerically shown to be robust across a wide range of macroscopic parameters of the full model. The wave speeds are thus found to depend strongly on the proliferation rate of the fittest population, resembling a winner-takes-all dynamics.

Species at the highest trophic levels are generally the least resilient to perturbations, yet the factors that determine their persistence remain poorly understood. We ask how the intrinsic self-regulation of a basal host governs the persistence of the trophic levels above it. To address this question, we use a minimal discrete-time model of a host-parasitoid-hyperparasitoid system motivated by the biology of trigonalid wasps, obligate hyperparasitoids that complete development only when their caterpillar host is also attacked by a primary parasitoid. We compare two forms of host self-regulation in the absence of natural enemies: scramble (Ricker) fecundity, appropriate when the parasitoid is the host's principal enemy and which recovers the classical model of Beddington and Hammond as a special case; and contest (Beverton-Holt) fecundity, appropriate when the host is held in check by many additional enemies. Using numerical continuation, Lyapunov exponents, and direct computation of the model's trajectories, we find that the form of host regulation is decisive for persistence at the top trophic level. Under scramble host dynamics, coexistence is extremely fragile with respect to both parameter values and initial conditions. Under contest host dynamics, coexistence is robust across broad parameter ranges and loses stability only gradually, through a supercritical Neimark--Sacker bifurcation and a subsequent torus-doubling route to chaos. Therefore, in at least some trophic networks, the resilience of the highest levels can be determined at the base of the network.

Background: Tuberculosis remains a persistent infectious disease whose control is complicated by latent infection, delayed treatment, incomplete recovery, reinfection, and continuing transmission from infectious individuals. Although treatment is central to tuberculosis management, it is frequently represented only as a transition parameter in mathematical models rather than as a separate epidemiological state. In this study, treatment was therefore incorporated explicitly as an independent compartment so that its influence on transmission, recovery, disease-induced mortality, and long-term disease persistence could be evaluated. Methods: A deterministic nonlinear compartmental model was formulated by dividing the total population into susceptible, exposed, actively infected, treated, and recovered classes. Reinfection of recovered individuals, progression from latent infection to active disease, movement of infectious individuals into treatment, treatment-associated recovery, natural mortality, and disease-induced mortality were included. Positivity and boundedness of the solutions were examined to establish biological validity. The basic reproduction number, R0, was derived through the next-generation matrix approach. Disease-free and endemic equilibria were determined, and their local and conditional global stability properties were investigated using Jacobian analysis, the Routh-Hurwitz criterion, center manifold theory, Lyapunov functions, and LaSalles invariance principle. Normalized sensitivity indices, Latin hypercube sampling, partial rank correlation coefficients, and numerical simulations were also applied. Results: The disease-free equilibrium was shown to be locally asymptotically stable when ,R0<1 whereas sustained transmission and a unique endemic equilibrium were associated with R0>1. Under the stated reduced-model assumptions, stability of the endemic equilibrium was established. Transmission-related parameters were identified as the strongest positive contributors to disease persistence. In contrast, treatment and recovery parameters were found to reduce the reproduction number and infectious burden. Numerical simulations indicated that stronger treatment implementation and reduced transmission opportunities produced substantial reductions in active tuberculosis cases. Conclusion: Treatment was shown to function as both a clinical pathway and an epidemiological control mechanism. The proposed framework may support the design of treatment-centered strategies for reducing tuberculosis prevalence and preventing long-term endemic persistence.

We introduce a minimal individual based model for population intelligence that retains three structures absent from standard single agent descriptions: a manifold valued intelligence state, endogenous population turnover, and heritable lineage dependent reproductive fitness. The validated S^1 intelligence diffusion was extended with neutral logistic birth death dynamics. Each individual occupies an angle on the circle S^1; noisy mean field attraction produces collective synchronization; logistic birth death dynamics regulate population size; and reproductive rates can depend on lineage and position relative to designated safe or unsafe directions. This nested construction is intentionally small enough to admit an analytic equilibrium benchmark while remaining capable of selection, invasion, and safety transitions. Population size returned toward carrying capacity from (N(0)=200,800,1600), with tail means (775.7,800.3,805.1) for K=800. Across K=100 to 1600, the intelligence equilibrium error scaled as K^{-0.513} , consistent with the canonical K^{-1/2} fluctuation rate, while the extinction rate was zero in 120 runs over the simulated horizon. The mean order parameter remained close to the original theoretical (r*=0.7242), as neutral demography predicts. Directional reproductive selection shifts a safety score to +0.7447 or -0.7448 while synchronization remains approximately 0.75, demonstrating that coordination is not equivalent to safety. These are controlled numerical results, not asymptotic proofs.

Ecological interactions play a crucial role in driving speciation and adaptive radiation. According to the adaptive dynamics theory of evolution in one-dimensional trait spaces, such evolutionary diversifications can occur when evolutionary branching points (convergence stable and evolutionarily unstable singular points) exist in the trait spaces. However, these evolutionary diversifications may occur across multiple traits at once, as exemplified by various within-species ecotypic diversifications involving multiple traits or genes. This paper attempt to extend this branching point condition to evolution in multi-dimensional trait spaces. Two criteria, termed branching possibility and branching inevitability, are developed to evaluate the likelihood of evolutionary branching in multi-dimensional trait spaces. Branching possibility guarantees the existence of evolutionary pathways that lead to branching (where deterministic parts of these pathways are described with the canonical equation of adaptive dynamics theory). Branching inevitability guarantees a steady progression of evolutionary branching (where the progression degree is described with a time-increasing function related to the local Lyapunov function, the potential function in an evolutionary potential game, and Fisher's fundamental theorem). Numerically simulated evolution in multi-dimensional trait spaces show that points having branching possibility or inevitability can induce evolutionary branching not only under asexual reproduction but also under sexual reproduction.

Biochemical Systems Theory (BST) represents nonlinear biochemical rate laws by local power-law approximations in logarithmic concentration coordinates. First-order coefficients are elasticities, whereas higher-order derivatives describe local log-synergism and its variation. Ordinary higher derivatives, however, are not tensorial under nonlinear reparameterizations and can mix biochemical response structure with coordinate artifacts. We formulate a covariant hierarchy, cBST1-cBST3, on a positive operating-point space equipped with a declared reference connection. cBST1 recovers classical elasticities in a flat logarithmic chart, cBST2 is the covariant Hessian of the log-response, and cBST3 is the symmetrized covariant derivative of cBST2. The framework is quantitatively evaluated using three representative rate laws from the curated yeast glycolysis model BIOMD0000000064: glucose transport, glucose phosphorylation, and phosphofructokinase. For 10,000 finite log-concentration perturbations at each of five radii, cBST2 reduced the cBST1 log-rate root-mean-square error by 96.6-98.6% at the largest tested radius, and cBST3 provided a further 68.4-98.1% reduction. The contracted cBST2 and cBST3 terms strongly predicted the corresponding lower-order truncation errors. Under the nonlinear transformation qi = sinh(ui), covariant contractions agreed across coordinates to within a 95th-percentile relative error of 1.2 x 10-13, whereas ordinary higher derivatives showed order-unity coordinate mismatches. Supplementary analytic tests recovered the expected second-, third-, and fourth-order truncation-error scaling. These results show that cBST1-cBST3 are not only coordinate-consistent descriptors but also practical diagnostics of where local power-law approximations require higher-order correction.

Ecological stoichiometry allows us to investigate how food quality affects food-web population dynamics. We consider periphyton, phytoplankton, and a shared Daphnia consumer in a closed system with separate benthic and pelagic phosphate pools. The producers have variable phosphorus-to-carbon ratios, while the consumer maintains a constant ratio. We construct the model, establish boundedness under stated conditions, characterize its grazer-free equilibrium families and their local spectra, derive conditions for positive consumer growth, and examine its numerical bifurcation structure as light-dependent carrying capacity varies. The bifurcation diagrams contain stable and unstable equilibrium branches together with a stable periodic family, revealing a more complex response to light enrichment than a single sequence of steady states and oscillations. A descending consumer equilibrium branch illustrates that greater light availability need not increase consumer biomass. We then introduce periodically varying carrying capacities and compare weak and strong seasonal forcing. Stronger forcing produces brief consumer rebounds in the high-light example, although the peaks diminish over the displayed interval. Together, the growth conditions and branch diagrams provide a framework for interpreting population dynamics in terms of food quantity and food quality. Independent numerical continuation identifies a Hopf crossing and a fold; the termination mechanism of the periodic branch remains unresolved.

Biological populations have been evolving in multi-dimensional trait spaces. Ecological interactions can induce the evolutionary diversification of a biological population into distinct multiple populations, which is called evolutionary branching in adaptive dynamics theory. If possible mutations are restricted to certain directions due to developmental, physiological, or physical constraints, the resulting evolutionary dynamics may be restricted to subspaces with fewer dimensions than the original trait spaces. In general, such constraints are also evolving slowly, due to rare and slight mutational deviations from the constraint subspaces. In this sense, constraints are not strict. However, the method for describing non-strict constraints that is globally consistent has yet to be developed. This study proposes a general and flexible method for constructing a non-strict constraint subspace from the corresponding strict constraint subspace. This method enables efficient analysis of adaptive evolution including evolutionary branching under non-strict constraints in the framework of adaptive dynamics theory. Applying this method to a resource competition model shows that the shape of the carrying capacity distribution can be treated as a non-strict constraint curve, and that the likelihood of evolutionary branching is affected not only by the shape of the constraint curve but also by its slow evolution.

Dengue transmission is shaped by multiple viral serotypes, temporary cross-immunity (TCI), antibody-dependent enhancement (ADE), and repeated exposure in endemic populations. Classical multi-strain models usually assume lifelong protection against reinfection with the same serotype. However, recent evidence suggests that homologous dengue reinfections, although rare, can occur. Their population-level consequences remain poorly understood. We extend a two-infection, two-strain dengue model with TCI and ADE-mediated transmission differences to include homologous reinfections. Homologous reinfection is represented by two exploratory parameters: relative susceptibility to reinfection with the same serotype and relative infectiousness during homologous reinfection. Using equilibrium analysis, bifurcation diagrams, simulations, and phase-space projections, we examine how these parameters affect dengue dynamics and interact with TCI duration and seasonal forcing under intermediate and long TCI durations, with and without seasonality. The extended model shows that qualitative dynamics characteristic of endemic dengue transmission are reproduced mainly when susceptibility to homologous reinfection is low, so that homologous reinfections remain rare but dynamically influential. Longer TCI broadens regions of complex oscillatory dynamics, while seasonality shifts the bifurcation structure and makes torus bifurcations a central route to complex behavior. Although backward bifurcation can occur when homologous susceptibility exceeds the biologically meaningful range, this result should be interpreted as a mathematical mechanism rather than a realistic dengue scenario. These results indicate that rare homologous reinfection pathways can influence long-term dengue dynamics when interacting with immune history, TCI, ADE-mediated transmission differences, and seasonal variation. Incorporating such pathways may improve understanding of recurrent outbreaks and irregular incidence patterns in highly exposed populations.